Lines

ImportantSummary

This tutorial explore how to handle spatial lines in R with terra package:

• read a spatial object with (terra::vect())
  
• calculate length of lines with (terra::perim())
• calculate distance to the nearest point with (sf::st_distance())

TipThe ecologist mind

Otter lives around rivers, but how far from the rivers were they observed? To answer this question, we will need to discover another type of vector: lines.

Setup

Follow the setup instructions if you haven’t followed the tutorial on points

If haven’t done it already, please follow the setup instructions. Let’s start with loading the required packages.

suppressPackageStartupMessages({
  library(mapview)
  library(here)
  library(sf)
  library(terra)
})

local

pt_otter <- vect(here("data", "gbif_otter_2021_mpl50km.gpkg"))

online

pt_otter <- vect(
  "https://github.com/FRBCesab/spatial-r/raw/main/data/gbif_otter_2021_mpl50km.gpkg"
)

Load lines from a geopackage file

Lines are made of multiple points. It is possible to create lines directly from coordinates but, in practice, it often comes from an existing spatial dataset. For our example, we will load rivers for the area of interest from IGN data BD CARTO.

Note that dataset has rough resolution (BD TOPO would be more suited for real analysis), but it’s perfect for our illustration and learning purpose.

You can load all vector dataset in R with the function sf::st_read(). It can read many different formats. Another function sf::read_sf() is also available and return a similar spatial object but with a tibble instead of a data.frame as attribute table.

terra

river <- vect(here("data", "BDCARTO-River_mpl50km.gpkg"))

sf

river_sf <- st_read(here("data", "BDCARTO-River_mpl50km.gpkg"))

Reading layer `BDCARTO-River_mpl50km' from data source 
  `/home/romain/GitHub/spatial-r/data/BDCARTO-River_mpl50km.gpkg' 
  using driver `GPKG'
Simple feature collection with 2110 features and 8 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 3.18983 ymin: 43.21278 xmax: 4.56447 ymax: 44.16754
Geodetic CRS:  WGS 84

online

river <- vect(
  "https://github.com/FRBCesab/spatial-r/raw/refs/heads/main/data/BDCARTO-River_mpl50km.gpkg"
)

NoteYour turn

1. How many different lines of river was loaded?
2. What is the coordinate reference system (CRS) of the loaded river data?
3. How can you get the name of the all the river stretches?

Click to see the answer

1. There was 2110 lines in the dataset. You can access it with dim(river), nrow(river) or just by typing river in the console.
   
2. The coordinates are in WGS 84 (EPSG 4326). You can access this information with crs(river, describe = TRUE) (or in sf with st_crs(river_sf)).
3. The column that stored the name of river stretches is toponyme. You could identify it with head(river) or names(river). There are 151 river stretches without names (table(is.na(river$toponyme))).

Calculate length of lines

terra

The function terra::perim() While calculating distance and length, be careful with projection systems. Some are not suited to calculate distance. The package terra recommends the calculation of distances in lat/long to get more accurate results (considering the geodesic distance, so accounting for Earth’s curvature).

# calculate the length of rivers (in km)
river$length_km <- perim(river) / 1000

sf

The function sf::st_length() calculate the length of lines. While calculating distance and length, be careful with projection systems. Some are not suited to calculate distance.It is recommended to use equidistant projections or use local projection system (if your study area is small).

# calculate the length of rivers (in km)
river_sf$length_km <- st_length(river_sf) / 1000

# see the distribution of river length
boxplot(river$length_km, ylab = "length (km)")

NoteYour turn

• Which is the longest river in our dataset?
  

Click to see the answer

# get the name of the longest river
river$toponyme[which.max(river$length_km)]

[1] "l'Hérault"

Map multiple layers

Let’s map the river and their length, as well as position of the otter observations.

interactive

You can combine multiple layers in mapview with a +.

mapview(river, zcol = "length_km") +
  mapview(pt_otter, col.regions = "red", color = NA)

terra

plot(river, y = "length_km")
plot(pt_otter, add = TRUE)

sf

plot(river_sf["length_km"], reset = FALSE, axes = TRUE)
plot(pt_otter, add = TRUE)

Warning

If you want to overlay multiple spatial object with base plot from sf, don’t forget to use the argument reset=FALSE.

Calculate the distance to the closest line

TipThe ecologist mind

At which distance from a river were the otter observed?

This question is addressed by the function terra::nearest().

Before comparing two spatial objects, it is recommended to plot the data (if not too big) and make sure the projection systems are the same and the extents match. In this case, do not use mapview (interactive map) because it will automatically project your data.

# make sure the spatial objects have the same projection
st_crs(river) == st_crs(pt_otter)

[1] TRUE

Find the closest river

terra

# Step 1: Find index of nearest line for each point
nearest_river <- nearest(pt_otter, river)

sf

In sf, this is a two-step process: (1) find the closest river sf::st_nearest_feature() and (2) calculate the distance between points and their closest river sf::st_distance().

# Step 1: Find index of nearest line for each point
nearest_sf <- st_nearest_feature(st_as_sf(pt_otter), river_sf)

The object nearest_river is a vector containing the index of the closest river for each otter observation.

# Step 2: Calculate distance only to the nearest line
nearest_distances <- st_distance(
  st_as_sf(pt_otter),
  river_sf[nearest_sf, ],
  by_element = TRUE
)

boxplot(nearest_distances, ylab = "distance to river (m)")

# should check out discrepancy
# plot(nearest_distances, nearest_river$distance)

NoteYour turn

1. Which are the rivers with most sights of otter? answer doesn’t work anymore
2. Make a map with the distance to river and visually check the coherence of the calculation.
3. (Reflexion) How can you explain the large distances?

Click to see the answer 1

table(river$toponyme[nearest_river]) |>
  sort(decreasing = TRUE) |>
  head(5)

Click to see the answer 2

# add the distance to river to the spatial object
pt_otter$dist_river <- nearest_distances

# make an interactive map
mapview(river) +
  mapview(pt_otter, z = "dist_river")

# or static map
plot(st_geometry(river))
plot(pt_otter[, "dist_river"], add = TRUE, pch = 16)